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We investigate band gaps, equilibrium structures and phase stabilities of several bulk polymorphs 
of wide gap oxide semiconductors ZnO, Ti 02 , Z 1 -O 2 and WO 3 . We are particularly concerned 
with assessing the performance of hybrid functionals built with the fraction of Hartree-Fock exact 
exchange obtained from the computed electronic dielectric constant of the material. We provide 
comparison with more standard density-functional theory and GW methods. We finally analyze 
the chemical reduction of Ti 02 into ThOa, involving a change in oxide stoichiometry. We show 
that the dielectric-dependent hybrid functional is generally good at reproducing both ground-state 
(lattice constants, phase stability sequences and reaction energies) and excited-state (photoemission 
gaps) properties within a single, fully ah initio framework. 
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I. INTRODUCTION 

Density-functional theory (DFT) has become the 
method of choice for calculating, understanding, and pos¬ 
sibly predicting properties of real materials from first- 
principles. The approach developed by Kohn and Shami 
(KS) provides a computational framework that allows 
to calculate ground-state related properties of realistic 
systems at affordable computational cost. For exam¬ 
ple, lattice constants of materials are often obtained in 
excellent agreement with experiment. In contrast, ex¬ 
cited states are usually poorly described, being beyond 
the scope of standard DFT-KS. A prototypical exempli¬ 
fication of this behavior is provided by the well-known 
failure of local and semilocal energy functionals in esti¬ 
mating the electronic band gap of semiconductors and 
insulators. In fact, the electronic gap obtained in the 
KS scheme (the KS gap) differs from true quasiparticle 
gap (measured e.g. in photoemission and inverse photoe¬ 
mission experiments) by the derivative discontinuity of 
the exchange-correlation (xc) energy functional at handj^ 
A rigorous description of electronic excitations requires 
working in the framework of many-body perturbation 
theory^ Within the latter, the GW approximation^ is to 
date the state-of-the-art method for computing quasipar¬ 
ticle band structures that can be meaningfully compared 
with experiments. However, GW calculations are much 


more computationally demanding than DFT-based ones, 
indeed limiting the size of the systems that can be inves¬ 
tigated at reasonable cost. Although methodological ad¬ 
vances have recently made it possible to treat remarkably 
complex surface/interface systems at the GW level^ 
DFT is still widely used for investigating electronic struc¬ 
ture features in bulk solids and surfaces, for which in most 
cases it is able to provide at least a qualitative picture. 
This can be justified on the basis of the formal similarity 
between the KS equation and the equation for the poles of 
the Green’s function of many-body perturbation theory, 
once the xc potential is identified with an approximation 
of the self-energy operator S. The development of hybrid 
functionals^ - — first proposed within the chemistry com¬ 
munity with the aim of increasing accuracy in molecular 
calculations and now routinely applied also to extended, 
condensed systems^ - — has allowed to gain even quan¬ 
titative accuracy in the computation of electronic band 
structures. 

Hybrid functionals are constructed by admixing a por¬ 
tion of Hartree-Fock exact exchange to a semilocal xc 
potential within the generalized-gradient approximation 
(GGA). Built upon a generalized KS scheme, the hy¬ 
brid functional approach provides a practical way to 
cure some of the known drawbacks of local and semilo¬ 
cal DFT xc functionals^! and constitutes an alterna¬ 
tive, cheaper electronic structure method with respect 
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to GW. This is especially true if localized basis sets 
are used for evaluating exact exchange, which allows to 
reduce the computational effort down to that of stan¬ 
dard calculations within the local-density approximation 
(LDA) or the GGA. Instead, using plane waves makes 
calculations about two order of magnitude more expen¬ 
sive, again limiting the size of the systems able to be 
addressed. Standard recipes for hybrid functionals pre¬ 
scribe a fixed portion a of exact exchange, indepen¬ 
dently of the system being investigated. This is the case 
of the full-range Perdew-Burke-Ernzerhof (PBE0)£ and 
screened-exchange Heyd-Scuseria-Ernzerhof (HSE06)^ 
hybrid functionals, in which a is set to 1/4. The latter 
has been identified as the optimal value for obtaining ac¬ 
curate atomization energies in molecular calculations 
While both PBEO and, above all, HSE06 have proved 
successful even when applied to solids,— it comes nat¬ 
urally to wonder if in this case a direct relationship 
between the exchange fraction a and some materials 
property can be established based on physical argu¬ 
ments. This possibility has been first investigated by 
Alkauskas et al. ^ who identified static electronic screen¬ 
ing as the essential ingredient determining the value 
of a in solids. This can be easily concluded by com¬ 
bining hybrid DFT and the static approximation to 
the GW self-energy, i.e. Hedin’s Coulomb-hole-plus- 
screened-exchange (COHSEX) approximation^ The ex¬ 
pression of the xc potential in the full-range hybrid func¬ 
tional scheme is 

v xc (r,r') = av^(r,r') + {l-a)v GGA {r) + v GGA {r) (1) 

where is the nonlocal Hartree-Fock exact exchange 
potential evaluated using KS orbitals, and w GGA and 
a GGA are the local GGA exchange and correlation po¬ 
tentials, respectively. The COHSEX self-energy can be 
partitioned as follows 

X(r, r') = S C oh ( r , r') + Esex(»", r ') ( 2 ) 

where the local Coulomb-hole (COH) and nonlocal 
screened-exchange (SEX) potentials are given by 

ScoH^r') = -\K r -r')[v(r,r') - W(r,r% 

N occ 

^SEx(r,r') = - y] <j>i{r)cj)i{r')W(r,r'). 

i= 1 

The screened Coulomb potential is defined through 

W(r,r') = J dr"e _1 (r, r")v(r", r') 

where v is the bare Coulomb potential and e” 1 is the 
inverse of the microscopic dielectric matrix, the evalua¬ 
tion of which represents one of the most computationally 
demanding steps in GW calculations. If microscopic con¬ 
tributions to dielectric screening are neglected and the in¬ 
verse dielectric matrix is approximated as the inverse of 


the macroscopic electronic dielectric constant e^ 1 , then 
the screened potential becomes W(r,r') « e^v(r,r'). 
Notice that the evaluation of the full microscopic dielec¬ 
tric matrix is now reduced to that of a ground-state, 
macroscopic, static quantity represented by the electronic 
dielectric constant e x . It is clear that such simplification 
implies considerable reduction in the computational ef¬ 
fort in comparison with GW. Finally, correspondence of 
spatially local and nonlocal contributions in Eqs. ID and 
ID leads to identifying 

a = 1/eqo, (3) 

where is the static electronic dielectric constant of 
the material— The corresponding hybrid xc functional 
is referred to as “dielectric-dependent” in the following. 

Since the paper of Alkauskas et al. several works have 
appeared in the literature adopting this paradigmatic 
scheme and evaluating Coo by means of different com¬ 
putational approaches—~« 7 Marques et aZ— computed 
Coo within the GGA Perdew-Burke-Ernzerhofi& approx¬ 
imation; alternatively, they suggested a relationship be¬ 
tween a and a density-dependent estimator of the band 
gap, which they used for defining a material-dependent 
hybrid functional. Skone et aZ— recently suggested a 
self-consistent scheme for obtaining a , by evaluating 
including local-field effects beyond the random-phase ap¬ 
proximation (RPA). Conesa— used instead the experi¬ 
mental dielectric constant. Other authors proposed ad¬ 
justing a so as to reproduce the experimental band gap 
of materials ) 18 ’ 2 ^ 30 aiming at testing how the functional 
performs in predicting different electronic features, such 
as band edge positioning at semiconductor interfaces and 
defect energy levels in solids. Concerning the last point, 
Chen and Pasquarello^l showed that good results are ob¬ 
tained if the parameter a is chosen in such a way that the 
modified hybrid functional matches the computed GW 
band gap of the intrinsic material. 

In all previous work in which a is derived through 
Eq. ©, attention has been focused on calculation of the 
electronic structure of semiconductors and insulators. In 
particular, the method has been shown to considerably 
improve over LDA, GGA and standard hybrid function¬ 
als concerning computation of the band gap for a wide 
class of materials. The authors of Refs. [2(1 [ID, and \2]\ 
also computed optimized lattice parameters for a set of 
solids, showing that an accuracy comparable to that of 
standard hybrid functionals is achieved. What has not 
been assessed so far is if the method is also capable to 
provide good results in terms of the energetics and ther¬ 
mochemistry of semiconducting materials. 

This is a crucial point to be addressed, both from a 
fundamental standpoint and in view of technological ap¬ 
plications. For example, in the field of catalysis and 
photocatalysis, metal oxide semiconductors are widely 
employed.— In a typical catalytic reaction, in which the 
semiconductor plays the role of the active catalyst, low- 
coordinated oxygen atoms at the semiconductor surface 
move into the reactant, leaving behind oxygen vacancies 
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in the original material, which are subsequently filled by 
molecular oxygen present in the atmosphere. Thus, the 
semiconductor undergoes a change in stoichiometry for 
which an accurate description of the energetics is de¬ 
terminant in order to properly predict its performance 
as a catalyst .— 3 The formation of a vacancy also possi¬ 
bly strongly affects its electronic properties. Accurately 
knowing the electronic band structure of the material is 
certainly paramount in the field of photocatalysis. For 
example, in semiconductor-based photocatalytic water 
splitting^ the band gap of the material has to exceed 
the redox potential of water (A-Ey^dox = 1.23eV), and 
in real devices should be no lower than 1.9 eV. Further¬ 
more, the valence band maximum and conduction band 
minimum of the semiconductor have to lie higher than 
the water reduction potential and lower than the water 
oxidation potential, respectively, in order to be able to 
reduce or oxidize water and produce molecular hydrogen 
or oxygen. Wide gap metal oxide semiconductors have 
been identified as promising candidate materials for pho¬ 
tocatalytic water splitting ; 32 ’ 35 ’ 36 as well as for catalysis, 
also due to their excellent optoelectronic properties. 

In this scenario, an assessment of the performance of 
the dielectric-dependent hybrid approach in predicting 
electronic and structural properties of oxide semiconduc¬ 
tors, as well as the energetics pertaining different kinds 
of transformations in such materials is certainly in order. 
We have selected a set of wide gap oxide semiconductors 
with current or potential technological impact, e.g. in 
the fields of catalysis, photocatalysis, and photovoltaics. 
As exemplified above, several issues concerning both the 
electronic structure and the energetics are to be tackled 
in these materials in order to properly conceive practical 
applications. To this purpose, it would be extremely use¬ 
ful to rely on a single theoretical framework in which all 
such properties can be accurately predicted from first- 
principles in a unified scheme. We have benchmarked 
the dielectric-dependent hybrid functional approach as 
a candidate method in this respect, allowing to reliably 
compute both band structures and total energies at rea¬ 
sonable computational cost. 

In the present work we systematically investigate band 
gaps and ground-state properties for the following set 
of oxide semiconductors in their bulk structures: ZnO, 
TiC> 2 , Zr 02 , WO 3 ; all are interesting materials for sev¬ 
eral applications in the energy field, as discussed above. 
Results are also reported for MgO, which is studied as a 
simple test material although being an insulator with a 
gap as large as 8 eV. DFT calculations have been per¬ 
formed using different local/semilocal and hybrid func¬ 
tionals, comparing plane waves and localized basis sets 
within the linear combination of atomic orbitals (LCAO) 
method. GW band gaps for selected polymorphs of the 
above oxides have also been computed for comparison. 
Then the relative stability of several polymorphs is in¬ 
vestigated at different levels of theory. Finally, the case 
of the chemical reduction of TiC >2 to TqOs, involving 
a change in oxide stoichiometry, is analyzed; it repre¬ 


sents a prototypical chemical process central to cataly¬ 
sis. The last two points require investigating, although 
from different perspectives, the common issue of the en¬ 
ergetics. For all of the above properties, we report a 
systematic comparison of at least one approach belong¬ 
ing to the classes of local/senrilocal, standard hybrid, and 
dielectric-dependent hybrid functionals. 


II. COMPUTATIONAL METHODS 
A. Plane waves 

Plane wave (PW) DFT calculations were performed 
using the QUANTUM ESPRESSO packaged A norm- 
conserving Troullier-Martins pseudopotential^ with 2s 
and 2 p electrons in the valence was used for oxy¬ 
gen, while for metals, norm-conserving Hartwigsen- 
Goedecker-Hutter pseudopotential s 39 ’ 40 explicitly includ¬ 
ing semicore electrons in the valence were employed. A 
careful treatment of core-valence partitioning in metal 
atoms has been proved to be determinant for obtaining 
reliable GW corrections^ especially in II-IV semicon¬ 
ductors, as analyzed in detail in Ref. [U It turns out 
that explicit treatment of semicore electrons of the metal 
atom in oxides is always needed, as we experienced in the 
case of ZnO (see discussion in Section IIII AD . Since our 
PW DFT calculations served as starting point for sub¬ 
sequent GW calculations, our choice of the pseudopo¬ 
tentials was guided by the above considerations. This 
means that (n — l)s and (n — 1 )p electrons of the metal 
atom were always explicitly treated in the valence, where 
n is the main quantum number of the outermost occu¬ 
pied electronic shell. Pseudopotentials were tested by 
checking that computed optimized lattice parameters and 
electronic band structures were in agreement with the 
reference literature, as well as that the former showed 
the expected trends at the LDA and GGA levels of the¬ 
ory when compared to experiment (see the Supplemental 
Material).— 

Including semicore electrons in the self-consistent cal¬ 
culation of charge density clearly requires using a large 
number of PW components in the expansion of the KS 
orbitals, as semicore electrons are tightly bound to the 
nucleus. This dramatically affects the resulting compu¬ 
tational cost, especially in hybrid functional calculations, 
in which evaluation of the nonlocal exact exchange using 
a PW basis set under periodic boundary conditions is 
particularly demanding. 

Brillouin zone (BZ) sampling is also critical for obtain¬ 
ing well-converged results. We chose Monkhorst-Pack^i 
(MP) k -point grids in such a way that estimated errors 
in computed band gaps and total energies were within 
20 meV. In hybrid calculations, evaluation of exact ex¬ 
change energy in PW basis requires the definition of an 
additional grid of g-points, which must be a symmetry¬ 
preserving subset of the fc-point grid. The choice of both 
k- and g-point grids determines the accuracy of the com- 
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puted quantities, while strongly affecting the computa¬ 
tional cost (see the Supplemental Material). 

Full structural optimizations were carried out only 
at the LDA and GGA levels, as optimizations with 
hybrid functionals are very expensive with QUANTUM 
ESPRESSO. A quasi-Newton Broyden-Fletcher-Goldfarb- 
Shanno (BFGS) scheme was adopted for the search of 
the energy minimum. Convergence thresholds on atomic 
forces and pressure were set to ICC 3 a.u. and 0.5 kbar, 
respectively. 


B. Localized basis sets 

The LCAO method, as implemented in the CRYSTAL09 
packaged was employed for DFT calculations using lo¬ 
calized basis sets. All-electron calculations were per¬ 
formed for MgO, ZnO and Ti02, while for Zr02 and WO 3 
the inner electrons of Zr and W atoms were described 
through effective core pseudopotentials (ECPs). In the 
CRYSTAL09 code, KS orbitals are expressed in terms of lo¬ 
calized atomic-like functions which are in turn expanded 
in a linear combination of Gaussian-type orbitals (GTOs) 
with fixed coefficients. Expansion coefficients and expo¬ 
nents of GTOs define a basis set for a given atom. 

In the present work, we adopted all-electron and va¬ 
lence basis sets which were previously tested in analo¬ 
gous solid state calculations. The following Gaussian- 
type all-electron basis sets were employed in our calcu¬ 
lations: 8-511 (dl) for Mg from Ref. [4(| pob-TZVP for 
Zn from Ref. E3, 8-4611(d41) for Ti from Ref. HJ For 
O the 8-411(dl) basis set of Ref. Hsj was used, except 
in the case of Z 1 -O 2 , in which the basis set from Ref. [50| 
was adopted, following the choice of a previous investi¬ 
gation on this material^ Small-core Hay-Wadt ECPs^ 
were employed for the heavy Zr and W atoms. For Zr 
only the 4 d and 5 sp electrons were explicitly treated, us¬ 
ing the 311(d31) valence electron basis set of Ref. l53l. For 
W, electrons belonging to shells 5 p, 5 d, 6 sp were included 
in the valence, and a modified Hay-Wadt double-^ basis 
set^ was employed to describe them. 

In the analysis of the chemical reduction of Ti02 to 
Ti 203 , for both Ti and O atoms, we adopted differ¬ 
ent basis sets than mentioned above, following previous 
Hartree-Fock and DFT+U studie s 55 ’ 56 on the same ma¬ 
terials. In molecular calculations we chose basis sets for 
oxygen-S and hydrogen^ such that experimental atom¬ 
ization energies for H 2 , O 2 and H 2 O molecules were well- 
reproduced at the PBEO level. 

As a high accuracy is needed in studying the energet¬ 
ics, BZ samplings were performed so as to ensure con¬ 
vergence within 1 meV on total energies for the different 
polymorphs of the materials under investigation. This 
corresponds to MP grids including a number of fc-points 
in the irreducible BZ at least equal to what we used in 
PW calculations (see the Supplemental Material). The 
lower cost of LCAO calculations made this choice afford¬ 
able. 


Different cutoff thresholds for the evaluation of 
Coulomb and exchange integrals written in terms of 
GTOs were set to their standard values in CRYSTAL09(^ 
The self-consistent field was considered converged when 
total energy difference between two subsequent cycles be¬ 
came lower than ICR 6 a.u. We kept these thresholds fixed 
throughout all the calculations for a more meaningful 
comparison of the computed total energies. 

Geometry optimizations were performed at both lo¬ 
cal/semilocal and hybrid functional level by fully opti¬ 
mizing both lattice parameters and atomic coordinates. 
A quasi-Newton algorithm with a BFGS Hessian up¬ 
dating scheme was adopted for the search of the en¬ 
ergy minimum. Geometry optimizations were termi¬ 
nated when the maximum component and root-mean- 
square (rms) of energy gradients became lower than 
0.00045 and 0.00030 a.u., respectively, and the maximum 
and rms atomic displacements were below 0.00180 and 
0.00120 a.u., respectively. 


C. GW calculations 

GW calculations were performed using the BERKE- 
LEYGW codej 59 ’ 60 The non self-consistent Go Wo scheme 
was adopted, in which both the electronic Green’s func¬ 
tion and the screened Coulomb potential were con¬ 
structed starting from DFT-KS eigenfunctions an d 
eigenvalues computed at the PBE level of theory 
(GoWo@PBE). Quasiparticle (QP) energies correspond¬ 
ing to band index n are evaluated on top of KS eigenval¬ 
ues using the first-order perturbation formula 

eST = + Mil S(£« p ) - |<® ■ (4) 

This implicit equation is solved by using Newton’s 
method to linearize the energy-dependent self-energy 
E (E) around the zero-order QP energy^ 

E° nk = E™ + (V£ fc s | E (E™) - v™* \<2) , 

finally yielding 


zt’QP ^ zt’O 
'Ilfc — 'Il/c 


cLE/dE , q _ „ K s\ 
1-dE/dE 1 nk nk> 


No further iteration was carried out on either the 
Green’s function or the self-energy. In fact, whether a 
self-consistent GW approach really improves computed 
band gaps over Go Wo is still matter of debated DFT cal¬ 
culations were performed at the experimental geometry 
using QUANTUM ESPRESSO. The dielectric function was 
evaluated at the RPA level, the frequency dependence 
being accounted for through the plasmon-pole approxi¬ 
mation of Hybertsen and Louie. 60 

In the BERKELEY GW code, all the quantities of interest 
are expressed in a PW basis set, evaluating them in recip¬ 
rocal space. Summations over a large number of empty 
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states are required for obtaining well-converged polariz¬ 
ability and Coulomb hole self-energy. This calls for ex¬ 
tensive convergence studies on various cutoff parameters 
entering calculations. The size of the static dielectric ma¬ 
trix in reciprocal space is controlled through the cutoff 
energy 2?®^. Additionally, the number of empty states in 
the expressions of the static RPA polarizability y, and of 
the Coulomb hole (CH) contribution to dynamical self¬ 
energy, Ech, should be determined. Polarizability typi¬ 
cally converges more rapidly than self-energy, for which 
we found that no less than 900 empty states were needed 
in our calculations (see the Supplemental Material). Fi¬ 
nally, the number of terms in the PW expansion of the 
xc potential, v xc , in Eq. (J4J is determined by the cut¬ 
off energy E™ t . Care should be taken when performing 
convergence tests, as not all cutoff parameters can be 
converged independently of each other. In fact, the size 
of the dielectric matrix in reciprocal space is not inde¬ 
pendent of the number of empty states included in Ech- 
In some previous work, this tricky behavior has not been 
considered with the due attention, leading to false con¬ 
vergence issues, as pointed out by Louie and coworkers£2 
for the particularly critical case of ZnO. 

In our calculations, careful convergence studies were 
carried out for all the materials, ensuring that accuracy 
within 50 meV on QP band gaps was reached. To this 
purpose, we closely followed the procedure described in 
Ref. [63l to which we refer for further details concerning 
the implementation of the GW scheme in BERKELEYGW, 
as well as for a description of the steps to be followed 
for assessing the convergence of calculated quantities. In 
Table S2 of the Supplemental Material cutoff parame¬ 
ters used in our calculations are reported. As they are 
not expected to depend strongly on the BZ sampling^ 
we estimated them using a minimal lxlxl grid. For 
the fully-converged GW calculations, BZ samplings were 
performed as in DFT PW calculations. 


D. DFT xc functionals 

In the first part of the work, we present our computed 
band gaps using different local/semilocal and hybrid xc 
functionals, comparing PW and LCAO approaches. We 
also compare band gaps computed using both dielectric- 
depended hybrid functional and GW approaches. 

In our local/semilocal calculations we used the 
Perdew-Zunger expression for the correlation energy^ 
for the LDA, and the Perdew-Burke-Ernzerhof (PBE) 
parametrizationi£ for the GGA. Hybrid functionals con¬ 
sidered are the full-range PBEO^ with exchange frac¬ 
tion a = 1/4 and the screened-exchange HSE06i£ with 
a = 1/4 and range of the screened exchange interaction 
set to its standard value uj = 0.11 a.u. -1 <££ We also tested 
the B3LYP hybrid functionalj&ill which combines a frac¬ 
tion a = 1/5 of Hartree-Fock exact exchange with the 
GGA Becke-Lee-Yang-Parr (BLYP) xc functional^ 

In hybrid calculations with QUANTUM ESPRESSO, we 


chose pseudopotentials generated at the closest GGA 
level, i.e. PBE for PBEO and HSE06, and BLYP for 
B3LYP. We are not aware, to date, of any publicly avail¬ 
able packages for pseudopotential generation with hy¬ 
brid functionals. Use of a Hartree-Fock instead of a 
GGA pseudopotential in a previous study on silicon^ 
showed that little effect should be observed on the com¬ 
puted band gap^ We obtained all pseudopotentials from 
the QUANTUM ESPRESSO pseudopotential library^ ex¬ 
cept for the Zn atom, for which the pseudopotential was 
converted from the cpmd format j 39 i 40 i 71 

We constructed dielectric-dependent hybrid function¬ 
als by defining the exchange fraction according to Eq. ©. 
Electronic dielectric constants were computed within the 
couple-perturbed Kohn-Sham (CPKS) method^— as 
implemented in the CRYSTAL09 code. This method ex¬ 
ploits first-order perturbation theory to calculate the 
polarizability, and fully includes local-field effects, as 
discussed in Ref. We evaluated both at the 
PBE and PBEO levels, thus defining the PBEOopbe and 
PBEOopbeo dielectric-dependent functionals. The inclu¬ 
sion of some degree of self-consistency in the calculation 
of a was also tested in certain cases. 


III. RESULTS AND DISCUSSION 

A. Band gap: local/semilocal, hybrid functionals 
and GW 

In this Section, the following materials are investi¬ 
gated: cubic (c) rocksalt MgO, wurtzite (wz) ZnO, 
anatase (a) Ti 02 , tetragonal (t) WO 3 , room-temperature 
7 -monoclinic ( 7 -m) WO 3 . We computed band gaps 
within DFT (with both PW and LCAO schemes) using 
different functionals, as well as within GW; results are 
reported and compared in Table [Q along with selected 
experimental results. For a more meaningful comparison 
of the different computational approaches, and due to 
the prohibitive cost of geometry optimizations using hy¬ 
brid functionals with QUANTUM ESPRESSO, calculations 
in this Section were performed at the experimental ge¬ 
ometry (see the Supplemental Material). 

We emphasize that comparison between theoretical 
calculations and experimental data is typically far from 
straightforward when it comes to spectroscopic proper¬ 
ties of semiconductors and insulators. From the exper¬ 
imental point of view, different techniques give access 
to different physical observables. In (direct and inverse) 
photoemission experiments, electrons are removed from 
or added to the solid (charged excitations), probing hence 
the energy difference between the N and the N± 1 parti¬ 
cles systems. From the theoretical point of view, such 
total energy differences are described by the poles of 
the one-particle Green’s function G , commonly named 
the QP energies. Thus, the measured QP energies or 
their difference (the fundamental gap) can be put in di¬ 
rect comparison with values obtained by solving the QP 
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eigenvalue problem , 6 - after devising an approximation to 
the self-energy operator E, for example within the GW 
scheme, or by approximating it with a DFT xc potential 
(in which case the QP energies effectively reduce to the 
KS eigenvalues). 

In contrast, optical measurements involve neutral ex¬ 
citations, which can be seen as creation of electron-hole 
pairs. The binding energy of electron-hole pairs is hence 
inherently included. The theoretical description of neu¬ 
tral excitations requires in principle working with the 
two-particle (electron-hole) Green’s function, and such 
an analysis lies beyond the scope of the present work. 
Nonetheless, by knowing or guessing the exciton binding 
energy, optical measurements can be useful for extrap¬ 
olating QP gaps when direct and inverse photoemission 
measurements are not available for the material at hand. 
This situation is not uncommon, as discussed in the fol¬ 
lowing for our investigated materials. 

As a second source of possible disagreement with ex¬ 
periment, we point out that our calculations, being per¬ 
formed at frozen ion positions, do not take into ac¬ 
count electron-phonon interaction, which generally leads 
to renormalization of the band gap, reducing it with re¬ 
spect to the zero-temperature QP gap * 81 ! 82 

Given the above considerations, and before discussing 
the results of our calculations, it is worthwhile to give an 
overview of the experimental situation. 


1. Experimental data 

The fundamental gap reported for MgO was extrapo¬ 
lated from thermoreflectance studies^ following estima¬ 
tion of the exciton binding energy. 

In the case of ZnO, different optical experiment s 76 ! 77 
agree on the value of our reported QP gap, which was 
extracted from optical experiments by detailed analysis 
of the excitonic levels. 

For anatase Ti 02 , the situation is more complicated, 
since, to our knowledge, no experimental data are avail¬ 
able for the fundamental gap. In Ref. [78j, analysis of 
temperature dependence of the absorption coefficient al¬ 
lowed to estimate the indirect optical gap to be 3.42 eV. 
For rutile Ti 02 , measurements of both the QP and the 
optical band gap are available. An optical band gap of 
~ 3.0 eV has been reported in several investigations ^ 8,88 
while combined photoemission and inverse photoemission 
experiments yielded for the direct fundamental band gap 
the values of 3.3 ± 0.5 eV (Ref. [13) and 3.6 ± 0.2 eV 
(Ref. [85} ). Assuming a similar relationship between fun¬ 
damental and optical gap holds for anatase, we could 
tentatively estimate its QP gap to be in the range 3.7 
4.0 eV, as also suggested in Refs. [H andl87l. 

The only available data for the pure phases of Z 1 -O 2 
are from optical absorption measurements , 79 from which 
a direct optical gap of 5.78 eV for the tetragonal phase 
was inferred. Photoemission studies yielded smaller mea¬ 
sured gaps of 5.68 eV (Ref.[Hj) and 5.5 eV (Ref.|89|); how¬ 


ever these values cannot be attributed to any of the three 
Zr 02 phases, as the probed samples were amorphous or 
polyscrystalline. 

For WO 3 , several ultraviolet direct/inverse photoemis¬ 
sion studie s 80 ! 90 ! 91 have reported a band gap in the range 
3.28-3.39 eV. Instead, absorption measurements*^ -- — 
have shown large data dispersion, with an optical band 
gap spanning the range 2.6-3.2 eV. 


2. DFT LDA, GGA and standard hybrids: comparison 
between PW and LG AO 

Results of our calculations confirm that, as expected, 
hybrid functionals provide substantial improvement over 
well-known LDA/GGA gap underestimation. In par¬ 
ticular, band gaps computed at the HSE06 level are 
in good agreement with experiment for a-TiC> 2 , t-ZrC >2 
and y-m-WOa, as also confirmed by previous calcula¬ 
tions on titaniu m 87 ! 99 ! 100 and tungsten^ oxides. For 
the same materials, we found that PBEO overestimates 
the band gap with respect to experiment, while B3LYP 
generally gives better results; this behavior is in agree¬ 
ment with previous investigations ! 51 ! 54 ’ 101- - *!^ Quantita¬ 
tive disagreement with the literature, when present, may 
be attributed to different choices of geometry (optimized 
vs. experimental, the latter being considered in the 
present part of the work). Independently of the adopted 
theoretical treatment, we obtained an indirect gap for a- 
Ti 02 and t-Zr 02 and a direct gap for the other materials. 
Both values are reported in Table [T| 

In the case of MgO and ZnO, we observe that hybrid 
functionals perform differently from what previously dis¬ 
cussed. MgO is an insulator with a wide band gap of 
nearly 8 eV; its computed gap, even using hybrid func¬ 
tionals, is smaller than the experimental value. In the 
case of ZnO, due to the strong gap underestimation at 
the LDA/GGA level, and contrary to what we found for 
the other semiconductors, PBEO seems to perform well, 
similarly to B3LYP, while HSE06 is not sufficient to cor¬ 
rect the PBE underestimation. Recent studies using hy¬ 
brid functional s 20 ’ 21 ! 104- -^ confirm the observed trends 
for MgO and ZnO. 

We finally comment on the comparison between re¬ 
sults obtained within PW and LCAO schemes. From 
Table [H it is found that discrepancies in computed band 
gaps range from ~ 0.05 to ~ 0.30 eV, most of them be¬ 
ing within 0.20 eV, depending on the material and the 
functional used. The closest correspondence is found at 
the LDA/GGA level, while with hybrid functionals differ¬ 
ences become more substantial. Different results are the 
consequence of different kinds of approximations charac¬ 
terizing the two computational schemes. Generation of 
reliable pseudopotentials on the one hand and construc¬ 
tion of accurate basis sets for solid state calculations on 
the other, are certainly two of the main critical points 
in this respect. For example, in the case of ZnO, using 
a Zn pseudopotential with only 3 d and 4s electrons in 
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TABLE I: Fundamental band gaps (in eV) of the various oxides computed within local/semilocal, standard hybrid, 
and Go Wo schemes, comparing PW and LCAO calculations. All calculations were performed at the experimental 
geometry. When the fundamental gap is indirect, direct band gaps are reported in parenthesis. Estimated errors are 

20 meV and 50 meV for the DFT and GW results, respectively. 



Type 

Method 

LDA 

PBE 

PBE0 

HSE06 

B3LYP 

G 0 IT 0 @PBE 

Expt. 

MgO 


PW 

4.63 

4.79 

7.18 

6.44 

7.04 

7.88 


rs 

LCAO 

4.86 

4.93 

7.35 

6.67 

7.07 


7.83^ 

ZnO 


PW 

0.78 

0.81 

3.18 

2.46 

2.90 

3.06 


wz 

LCAO 

0.92 

1.02 

3.35 

2.75 

3.05 


3.44^ 

Ti0 2 


PW 

2.03 (2.16) 

2.13 (2.26) 

4.40 (4.47) 

3.67 (3.74) 

3.94 (4.03) 

3.73 (3.81) 

3.42^ 


LCAO 

2.09 (2.18) 

2.18 (2.28) 

4.23 (4.33) 

3.59 (3.70) 

3.79 (3.89) 


Zr0 2 


PW 

4.01 (4.18) 

4.12 (4.29) 

6.52 (6.69) 

5.78 (5.95) 

6.09 (6.25) 

5.87 (6.06) 

5.78^ 


LCAO 

3.94 (4.09) 

4.00 (4.16) 

6.33 (6.45) 

5.61 (5.75) 

5.90 (6.02) 


W0 3 


PW 

1.84 

1.94 

3.90 

3.16 

3.53 

3.34 


7 -m 

LCAO 

1.86 

1.91 

3.74 

3.09 

3.39 


3.38^ 


a Thermoreflectance at 85 K, Ref. [75 

b Transmission spectroscopy at 4.2 K, Ref. [76; photoluminiscence at 2 K, Ref. |77| 
0 Transmission spectroscopy, Ref. lid 

d Reflectance spectroscopy (VUV) in Y-doped Zr0 2 , Ref. l79l 
6 Photoemission (UPS-IPES), Ref. l80l 


valence, and adopting a non optimal basis set for Zn in 
LCAO calculations, resulted in computed gaps differing 
by as much as ~ 0.5 eV within LDA/GGA and ~ 1 eV 
at the hybrid functional level. Adoption of a better pseu¬ 
dopotential and basis set indeed partially corrected this 
discrepancy. Inclusion of semicore electrons in Zn has 
also been found to be crucial for obtaining reliable Go Wo 
corrections over the PBE gap^™ similarly to what was 
reported by Gori et alr^ In conclusion, we believe that, 
given the above discussion, the agreement between our 
computed PW and LCAO band gaps should be deemed 
satisfactory. 


3. Many-body perturbation theory: Go Wo 

Table Q] reports Go Wo band gaps computed on top of 
PBE calculations. We notice that generally good agree¬ 
ment is obtained with experiment for all the materials. 
Our results are in line with previous investigations, al¬ 
though numerical values may differ as a consequence of 
the different computational setup adopted, e.g. the start¬ 
ing DFT approximation to QP energies and wavefunc- 
tions, and the way the frequency dependence of the di¬ 
electric function is accounted for 

For MgO, our computed gap of 7.88 eV is very close 
to experiment and to the results of the Go Wo investiga¬ 
tion of Fuchs et al. adoping a HSE03 starting point 
Calculations starting from DFT-PBE using the projected 
augmented wave (PAW) method gave smaller QP gaps of 
7.25 eV (Ref. HU) and 7.41 eV (Ref. ITIol). 

For TiC> 2 , our obtained indirect QP gap of 3.73 eV 


is identical to the full-frequency result of Land- 
mann et al., 87 and in close agreement with other 
investigations 

Quasiparticle band structure calculations of ZrC >2 are 
scarce. The only Go Wo studies we are awar e of reported 
a fundamental band gap of 6.40 eV (Ref. fill and 5.56 eV 
(Ref. m for tetragonal zirconia. Our result of 6.06 eV 
for the direct gap is in between these two values, and 
compatible with the measured optical gap of 5.87 eV. 

As for 7 -UL-WO 3 , recent studies of Galli and cowork¬ 
ers y ielde d QP gaps of 3.26 eV (Ref. Ill 6 t ) and 3.30 eV 
fRef. Ill7h . which are very close to our result, and in ex¬ 
cellent agreement with the measured photoemission gap. 

Finally, ZnO stands out as a particularly critical case, 
for which considerable disagreement is found in the the¬ 
oretical literature concerning its QP properties. Sev¬ 
eral studies have yielded QP band gaps ranging from 
~ 2.1 eV (Refs. ( 6 l] and 1 1091 ) to 3.4 eV obtained by 
Louie and coworkers^ within Go IFo@LDA. These au¬ 
thors suggested that smaller computed gaps may be the 
consequence of false convergence behavior of the dielec¬ 
tric function and self-energy. Stankovski et studied 
the effect of different plasmon-pole models on the band 
gap, while Friedrich et air — carried out Go Wo on top 
of all-electron DFT calculations to investigate the effect 
of the pseudopotential approximation. At present there 
is no general consensus onto which is the most correct 
approach to describe the QP band structure of this ma¬ 
terial. Our calculations gave a gap significantly smaller 
than the experiment, which is in line with several pre¬ 
vious investigations within the Go W'o@LDA/GGA ap¬ 
proach. 












TABLE II: Band gap energy (in eV) computed using dielectric-dependent hybrid functionals defined through the 
exchange fractions reported in Table Hill Comparison with standard PBEO and Go Wo is provided. Only the direct 
gap is reported for Zr02, since the measured gap is extrapolated via a model for a direct gap. MAE and MARE are 
the mean absolute and mean absolute relative error, respectively, with respect to the computed Go Wo band gap. 
All DFT calculations were performed within the LCAO scheme at the experimental geometry. Estimated errors are 

20 meV and 50 meV for the DFT and GW results, respectively. 



Type 

PBEOqipbe 

PBEOopbeo 

PBEO 

Go lFo@PBE 

Expt. 

MgO 

rs 

8.06 

8.33 

7.38 

7.88 

7.83 

ZnO 

wz 

3.18 

3.94 

3.35 

3.06 

3.44 

Ti0 2 

a 

3.38 

3.72 

4.24 

3.73 

3.42 

Zr0 2 

t 

5.81 

6.12 

6.45 

6.06 

5.78 

WO 3 

7 -m 

3.23 

3.50 

3.74 

3.34 

3.38 

MAE (eV) 


0.20 

0.31 

0.42 



MARE (%) 


4.6 

8.1 

9.6 




TABLE III: Electronic dielectric constant 
(estimated error within 0.02) and corresponding 
exchange fraction a (%) evaluated at the PBE and 
PBEO levels. Calculations were performed within the 
LCAO scheme at the experimental geometry. Measured 
dielectric constants are provided for comparison. 



Type 

PBE 

PBEO 

Expt.^ 

€00 

a 

^OO 

a 

€00 

MgO 

rs 

3.10 

32.2 

2.87 

34.9 

2.96 

ZnO 

wz 

4.31 

23.2 

3.24 

30.8 

3.74 

Ti0 2 

a 

6.52 

15.3 

5.16 

19.4 

5.62 

Zr0 2 

t 

5.45 

18.3 

4.62 

21.6 

4.9 

WOs 

7 -m 

5.43 

18.4 

4.56 

21.9 

4.81 


a Measured static electronic dielectric constants from the 
following references: MgO, Ref. I 9 S Z nO. Ref. [ 9 3: Ti02, 
Ref. [ 93 J Zr0 2 , Ref. El, WO 3 , ReflU 


4- DFT dielectric-dependent hybrid functionals 

We turn to study the performance of different kinds 
of dielectric-dependent hybrid functionals. The fraction 
of exact exchange is defined through a = l/e^; the 
static electronic dielectric tensor was computed within 
both PBE and PBEO; £oo is defined as the average of the 
diagonal elements of the diagonalized dielectric tensor 
written in real space. 

Table EH reports the computed dielectric constants 
along with the corresponding exchange fractions apsE 
and apBEO; defining the modified PBEO functionals, 
PBEOapBE and PBEOapBEo- We notice that, as ex¬ 
pected, PBE and PBEO systematically overestimates and 
underestimates, respectively, the dielectric constant, as a 
consequence of the opposite tendency concerning estima¬ 
tion of the band gap. 

Table E! shows that the dielectric-dependent hybrid 
functional approach remarkably improves the computed 
gap over standard fixed-a PBEO. Agreement with exper¬ 
iment and Go Wq calculations is overall satisfactory with 


both modified hybrids, with the exception of ZnO, for 
which standard PBEO already performs well. Since reli¬ 
able experimental photoemission gaps are not available 
for all the materials, mean errors in Table El are defined 
with respect to the computed Go Wo gaps. Error anal¬ 
ysis shows that the PBEOapsE functional achieves the 
closest agreement with Go Wq for the class of materials 
under investigation. We thus took it as our preferred 
starting dielectric-dependent hybrid approach in the rest 
of the work. 


B. Oxide polymorphs: band gaps, crystal 
structures and phase stability 

We now turn to investigation of several other struc¬ 
tural modifications of our set of oxide semiconductors. 
This includes analysis of the their relative phase stabil¬ 
ity, electronic and dielectric properties, and equilibrium 
crystal structure, by using different theoretical methods. 

Phase stability is evaluated on the basis of ground- 
state total energy calculations, thus not taking into ac¬ 
count thermal and entropic contributions. In particular, 
the thermal contribution mainly comes from lattice vi¬ 
brations (phonons), and may be relevant in temperature- 
induced structural phase transitions. Here we argue that 
an analysis of the only zero-temperature electronic con¬ 
tribution is still able to provide a qualitative picture on 
the stability issue. In fact, although in this case quantita¬ 
tive comparison with thermochemical data may be of lim¬ 
ited significance, it can be interesting to assess whether 
theory is able to correctly reproduce the experimentally 
found phase stability sequence (at ambient pressure) for 
a given material. 

Since typically tiny energy differences are involved in 
transformations between phases with similar stabilities, 
properly reproducing the stability order is an extremely 
challenging task, pushing to the limits of accuracy of 
DFT methods. Different functionals may predict dif¬ 
ferent stability sequences, as a consequence of changes 
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TABLE IV: Electronic dielectric constant 
(estimated error within 0.02) and corresponding 
exchange fraction a = 1/cqo (%) evaluated within PBE 
and dielectric-dependent PBEO (PBEOopbe and 
PBE0apB E ) for various polymorphs of the studied 
oxides. Calculations were performed within the LCAO 
scheme at the optimized geometry. 



Type 

PBE 

PBEOapBE 

PBE0a™ B 


Coo 

a 

Coo 

a 

Coo 

a 

MgO 

rs 

3.13 

31.9 

2.80 

35.7 

2.78 

36.0 

ZnO 

wz 

4.30 

23.2 

3.28 

30.8 

3.14 

31.8 

zb 

4.67 

21.4 

3.37 

29.6 

3.19 

31.3 


r 

7.98 

12.5 

6.76 

14.8 

6.59 

15.2 

Ti0 2 

a 

6.58 

15.2 

5.59 

17.9 

5.45 

18.3 


b 

6.95 

14.4 

5.99 

16.7 

5.87 

17.0 


m 

5.29 

18.9 

4.62 

21.6 

4.55 

22.0 

Zr0 2 

t 

5.56 

18.0 

4.89 

20.5 

4.82 

20.8 


c 

5.86 

17.1 

5.09 

19.7 

4.99 

20.0 


e-m 

5.58 

17.9 

4.74 

21.1 

4.65 

21.5 


tr 

5.54 

18.1 

4.68 

21.4 

4.56 

21.9 

W0 3 

7-m 

5.55 

18.0 

4.68 

21.4 

4.57 

21.9 

or 

5.43 

18.4 

4.55 

22.0 

4.45 

22.5 


t 

6.47 

15.4 

5.33 

18.7 

5.17 

19.4 


c 

10.09 

9.90 

8.27 

12.1 

7.98 

12.5 


in relative total energies of the order of some meV per 
atom, going from one level of theory to another. 

In the following, we test the performance of the PBE, 
PBEO and PBEOojpbe functionals in this respect. Equi¬ 
librium geometries were also investigated, as it is ex¬ 
pected that a good xc functional is able to yield opti¬ 
mized lattice parameters within few percent from exper¬ 
iment. Moreover, accurate geometry optimizations are 
needed for obtaining meaningful total energies to be com¬ 
pared for analysis of phase stability. For this reason, in 
the present Section all the reported quantities (includ¬ 
ing dielectric constants, which in turn affect the amount 
of exact exchange entering the definition of the modified 
PBEO functional) were computed at the optimized geom¬ 
etry of the relevant xc functional. Since the CRYSTAL09 
code is very performant when it comes to structural opti¬ 
mizations with hybrid functionals, we used it throughout 
all the calculations presented in the remaining part of the 
work. 

In Table |W] the electronic dielectric constant and the 
corresponding exchange fraction calculated at different 
levels of theory are reported for the phases under investi¬ 
gation. A partial self-consistency on a = 1/cqo was car¬ 
ried out by evaluating within PBEOopbe, thus defin¬ 
ing the PBE0apg E functional. Comparison of Tables Mil 
and m shows the effect of geometry optimization on the 
computed dielectric constants for the commonly studied 


TABLE V: Differences in total energy per formula unit 
(in meV, estimated error 1 meV) with respect to a 
chosen phase (redefining the zero of the energy for each 
functional individually), for various polymorphs of the 
studied oxides, computed at different levels of theory. 



Type 

PBE 

PBEO 

PBEOc^pbe 

PBE 0 o'pg E 

ZnO 

WZ 

0 

0 

0 

0 

zb 

49 

54 

102 

80 


r 

20 

62 

18 

22 

Ti0 2 

a 

0 

0 

0 

0 


b 

-5 

32 

1 

2 


m 

-20 

17 

33 

25 

Zr0 2 

t 

0 

0 

0 

0 


c 

74 

45 

40 

36 


e-m 

5 

1 

-2 

-8 


tr 

0 

0 

0 

0 

wo 3 

7 -m 

or 

3 

6 

6 

10 

2 

36 

6 

1 


t 

30 

70 

-153 

-139 


c 

201 

301 

-403 

-462 


polymorphs at the PBE level. 

1. Polymorphs structures and relative stability 

Apart from wurtzite (wz), which is the most thermo¬ 
dynamically stable phase of ZnO at ambient conditions, 
the zinc-blende (zb) modification can be stabilized upon 
growth in specific conditions =82 The rocksalt structure 
exists only at high pressures, and it is not addressed here. 
From Table m it is seen that the better stability of the 
wurtzite phase is correctly predicted at all levels of the¬ 
ory, with the partially self-consistent PBEO confirming 
the result of the PBEOopbe- These findings are in agree¬ 
ment with previous calculations using LDA, GGA and 
screened-exchange hybrid functionals^ 21 

The three naturally-occurring polymorphs of Ti02 are 
rutile (r), anatase (a) and brookite (b). Their rela¬ 
tive thermodynamic stability critically depends on crys¬ 
tal size. At ambient conditions, the thermodynami¬ 
cally most stable macrocrystalline (crystal size exceed¬ 
ing 35 nm) phase of Ti02 is rutile. From calorimet¬ 
ric measurements of the phase transformation of anatase 
and brookite to rutile it is suggested that the phase sta¬ 
bility sequence is rutile > brookite > anatase^ 22 It is 
well-known that GGA is not able to reproduce the cor¬ 
rect sequence, predicting anatase to be more stable than 
rutile*i2^ Our results confirm such tendency and suggest 
that, independently of the amount of Hartree-Fock ex¬ 
change used, hybrid functionals are similarly not able to 
capture the correct sequence, in agreement with the pre¬ 
vious investigation of Labat et al > — Recent studie s 125 ’ 126 
demonstrated that dispersion van der Waals interactions 
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FIG. 1: Variation of total energy as a function of the exchange fraction a entering the definition of the 
dielectric-dependent hybrid functional. Total energies per unit MO x (M = Zn, Ti, Zr, W) are given with reference 
to that of a chosen polymorph for each material: (a) wurtzite for ZnO, (b) anatase for TiC> 2 , (c) tetragonal for 
Zr 02 , (d) 7 -monoclinic for WO3. In the insets convergence of the electronic dielectric constant Coo is analyzed. 


have to be taken into account in order to correctly repro¬ 
duce the experimental sequence. 

Zirconia features three different structural modifica¬ 
tions at ambient pressure; the baddeleyite structure, with 
a monoclinic (m) unit cell, is stable at room temper¬ 
ature, and transforms into a tetragonal (t) phase at 
1480 K, which is stable until 2650 K when it is converted 
into a cubic (c) flourite phase. In fact, the tetragonal 
and monoclinic phases can be viewed as distorted cu¬ 
bic structures^ From Table 0 it is inferred that while 
PBE predicts the correct stability sequence as a function 
of increasing temperature, PBEO inverts the monoclinic 
and tetragonal phases, in fact predicting the latter to be 
more stable than the former, in agreement with previous 
studies^ but opposite to experimental evidence. Varying 
the exchange fraction within the dielectric-dependent hy¬ 
brid scheme did not yield qualitatively different results, 
similarly to what was found for Ti 02 - 

Finally, the case of WO 3 deserves particular attention, 
since its phase diagram at ambient pressure is exception¬ 
ally rich, exhibiting five different polymorphs transform¬ 
ing into each other by variation of the temperature. The 


simple cubic (c) structure, only recently successfully sta¬ 
bilized at ambient conditions ; 12 - is constituted by a lat¬ 
tice of regular corner-sharing WC >6 octahedra, forming 
a strongly ionic compound. Several structures are suc¬ 
cessively obtained by decreasing the temperature from 
the melting point at 1700 K: tetragonal (t, stable above 
740 °C), orthorombic (or, from 330 °C to 740 °C), room- 
temperature 7 -monoclinic ( 7 -m, from 17 °C to 330 °C), 
triclinic (tr, from — 50 °C to 17 °C), and low-temperature 
e-monoclinic (e-m, from —140 °C to — 50°C)<^ 

Octahedra become more and more distorted going from 
the tetragonal to the low-temperature monoclinic phase. 
Previous work s 54 ' 117 ’ 128 investigated how such distortions 
affect the electronic band structure, as a consequence of 
energy-lowering rearrangement of the WOg units. The 
recent study by Ping and Galliiil suggested that tilting 
of the W0 6 octahedra along different crystalline axis is 
the main responsible for band gap opening going from the 
cubic and tetragonal phases (where no tilting is present) 
to the other, less symmetric structures. This is apparent 
from Table fVTl in which the reported band gap increases 
by as much as ~ 1 eV or more passing from the tetragonal 
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TABLE VI: Band gap energy (in eV, estimated error 20 meV) computed within PBE and dielectric-dependent 
PBEO for various polymorphs of the studied oxides. Fundamental and direct gaps are reported outside and inside 
parenthesis, respectively. Calculations were performed within the LCAO scheme at the optimized geometry. 



Type 

PBE 

PBEO 

PBEOapBE 

PBE0apg E 

Expt.^ 

MgO 

rs 

4.63 

7.38 

8.24 

8.67 

7.83" 

ZnO 

wz 

1.07 

3.50 

3.32 

4.07 

3.44 

zb 

0.91 

3.29 

2.93 

3.76 



r 

1.71 

3.90 

2.76 

2.96 

3.3 

Ti0 2 

a 

2.14 (2.23) 

4.36 (4.36) 

3.42 (3.46) 

3.66 (3.70) 

3.42 


b 

2.36 

4.46 

3.53 

3.71 



m 

3.38 (3.78) 

5.67 (6.10) 

5.09 (5.59) 

5.34 (5.77) 

5.83 

Zr02 

t 

3.80 (3.86) 

6.07 (6.20) 

5.43 (5.53) 

5.66 (5.76) 

5.78 


c 

3.12 (3.63) 

5.43 (6.04) 

4.67 (5.25) 

4.91 (5.50) 

6.1 


e-m 

1.80 

3.85 

3.28 

3.52 



tr 

1.67 

3.68 

3.06 

3.36 


W0 3 

7-m 

1.62 

3.67 

3.07 

3.34 

2.6 

or 

1.39 

3.34 

2.81 

3.09 



t 

0.53 

2.23 

1.53 

1.76 

1.8 


c 

0.55 (1.58) 

2.19 (3.21) 

1.15 (2.17) 

1.29 (2.31) 



a Measured photoemission gap of r-Ti02 from Ref. l84i; direct optical gap of m-, t-, c-ZrC >2 from Ref. |79>: optical gap of t- and 7 -m-W 03 
from Ref.[ll7. For the remaining references on measured band gaps see Table [Tj 


and cubic structures to the more distorted ones. A simi¬ 
lar relationship between structural and electronic proper¬ 
ties was found in a recent DFT study on doped W 03 r^ 
in which intercaled atoms drove the structural deforma¬ 
tion responsible for the large band gap reduction. 

As expected, the electronic dielectric constant follows 
the opposite trend with respect to the band gap, being 
larger in the cubic and tetragonal phases and decreasing 
with increasing band gap, as shown in Table IIVI This in 
turn affects the corresponding exchange fractions, whose 
strong dependence upon the phase considered is at the 
origin of the wrong prediction of the stability order of 
WO 3 within the dielectric-dependent hybrid scheme. In 
fact, Table El shows that, while using standard, fixed-a 
PBEO functional the cubic and tetragonal phases were 
predicted to be the least stable, in agreement with the 
experimental order as a function of temperature, they 
turn to be the most stable ones when a is obtained from 
the computed dielectric constant. This behavior is not 
observed in the other materials, where the electronic and 
dielectric properties vary more smoothly upon the phase 
transition, or, equivalently, the electronic structure cor¬ 
relates more weakly to lattice distortions. In those cases, 
the dielectric-dependent hybrid approach proved to qual¬ 
itatively reproduce the results of standard PBEO. This 
is clearly seen in Fig. [TJ where relative total energy is 
plotted as a function of the exchange fraction. Non¬ 
crossing lines on the positive energy half-plane indicate 
that no phase stability inversion occurs going from PBEO 
to PBEOapBE, contrary to what happens in the case of 
W0 3 . 


2. Band gap dependence on polymorphic structure 

The interesting case of WO 3 shows that the electronic 
structure can be significantly affected by the polymor¬ 
phic structure. Having assessed the performance of the 
dielectric-dependent hybrid method in band gaps com¬ 
putation of some well-characterized oxide polymorphs, it 
is worthwhile to investigate the electronic structure of 
other less common phases. 

Table IVD collects computed band gaps for all the phases 
studied so far. Here, all calculations were performed at 
the optimized geometries. Due to the experimental diffi¬ 
culty in obtaining samples of a well-defined phase, mea¬ 
sured gaps are not available for all the structures or are 
very disperse, as in the case of brookite Ti 02 - 13 ° Hence, 
comparison with previous hybrid and GW results pro¬ 
vides a good benchmark for our calculations. 

In general, the PBEOapBE confirms itself as an ex¬ 
cellent method for computing accurate electronic struc¬ 
tures. Notice however that self-consistency on the di¬ 
electric constant does not always improve the computed 
gap, as happens in the case of MgO and ZnO. Analysis 
of Table IIVI and of the insets in Fig. [T] shows that con¬ 
vergence on the dielectric constant is rapidly achieved, as 
also reported in Ref. US Due to the considerable cost of 
geometry optimizations for the largest cells, we limited 
ourself to a single self-consistency step on a, being confi¬ 
dent that PBECtapgg results are not far from that of the 
fully self-consistent PBEO. 

In the following, we briefly compare our results with 
the existing theoretical and, when available, experimen- 
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TABLE VII: Optimized cell parameters for various 
polymorphs of the studied oxide semiconductors, 
computed at different levels of theory. Calculations were 
performed within the LCAO scheme using CRYSTAL09. 



Type 

Parameter 

PBE 

PBEO 

PBE 0 apg E 

Expt.^ 

MgO 

rs 

a (A) 

4.258 

4.207 

4.185 

4.212 

ZnO 

wz 

a (A) 
c (A) 

3.292 

5.119 

3.264 

5.077 

3.257 

5.071 

3.250 

5.207 

Ti0 2 

a 

a (A) 
c (A) 

3.821 

9.672 

3.765 

9.656 

3.789 

9.631 

3.781 

9.515 

Zr0 2 


a (A) 

3.640 

3.608 

3.613 

3.571 


c (A) 

5.288 

5.209 

5.224 

5.182 



a (A) 

7.444 

7.334 

7.348 

7.306 

W0 3 

7 -m 

b (A) 
c (A) 

7.672 

7.885 

7.605 

7.801 

7.616 

7.817 

7.540 

7.692 



p n 

90.65 

90.61 

90.62 

90.88 


a See the Supplemental Material. 


tal literature. 

For ZnO, previous hybrid calculations^ predicted the 
band gap in the zinc-blende phase to be ~ 0.2 eV smaller 
than in wurtzite, which is also supported by experimental 
evidencefi^i confirming our PBEOapeE results. 

For Ti 02 we obtained a direct gap for rutile which 
is smaller by ~ 0.7 eV than in the anatase phase. 
Our computed gap for rutile is smaller than previous 
Go IEo@LDA/GGA results spanning the range from 3.3 
to 3.6 e V t 86 i 87 i 112 i 132 but close to the value of 2.85 eV ob¬ 
tained within Gq ITo@LDA-|-Uri^ and compatible with 
the measured QP gap of 3.3 ± 0.5 eV.— For the brookite 
phase, a recent Go W 0 @PBE investigation^ yielded a 
direct band gap of 3.86 eV, larger than in rutile and 
anatase, and in good agreement with our PBEOapgg cal¬ 
culations. 

For Zr 02 , our predicted band gap ordering for the 
three phases is in line with previous investigations us¬ 
ing hybrid functionals^ and GIF ) 115 : 133 and numeri¬ 
cal val ues a gree very well with the GIFo@LDA results 
of Ref. Ill5l reporting an indirect band gap of 5.34 eV, 
5.92 eV and 4.97 eV for the monoclinic, tetragonal and 
cubic structures, respectively. 

Only recently GW calculations were performed for all 
the ambient pressure polymorphs of WO 31 117 Computed 
band gaps within Go IFo@LDA are in good agreement 
with our PBECtapgp results, and confirm our obtained 
band gap ordering. 

3. Equilibrium geometries 

Finally, in Table IVIII we report the optimized cell 
parameters computed within PBE, standard PBEO 


and dielectric-dependent PBEO for some selected poly¬ 
morphs. Full optimization of both cell parameters and 
atomic positions was carried out. Results indicate that 
dielectric-dependent PBEO generally performs as well as 
standard PBEO when it comes to structural properties, 
the error with respect to experimental lattice constants 
being in most cases within ~ 2%. 


C. Stoichiometry transformations: the chemical 
reduction of Ti0 2 

We finally assess the performance of the dielectric- 
dependent hybrid method in calculating reaction energies 
for the chemical reduction of TiC >2 to Ti 2 03 . Differences 
in total energy involved in chemical reactions are typi¬ 
cally of the order of few eV, compared to few meV or 
tens of meV, observed in structural phase transitions of 
the previous section. This allows us to meaningfully ad¬ 
dress the issue of whether the method is able or not to 
quantitatively improve over standard hybrid functionals 
in terms of computed total energies. 

We are concerned here with the two reduction reactions 
of Ti0 2 : 

2Ti0 2 + H 2 —» Ti 2 0 3 + H 2 0 (5) 

2 Ti0 2 —> Ti 2 03 + — 0 2 . (6) 

TABLE VIII: Relative total energies E tot (in meV, 
estimated error 1 meV) per unit Ti 2 0 3 , band gap E g (in 
eV, estimated error 20 meV) and lattice constants a and 
c (in A), computed at the PBE0apg EO level. AFM 1 and 
AFM2 labels two different antiferromagnetic (AFM) 
configurations; FM refers to the ferromagnetic one. 
a = 12.8% was used in the dielectric-dependent PBEO; 
this value was obtained from the electronic dielectric 
constant 6 ^ = 7.81 computed within PBEOapsEo in the 
AFM1 configuration. Total energies are reported 
relative to that of the lowest-energy AFM1 phase. For 
spin configurations the notation (Tii, Ti 2 , Ti 3 , TL), 
with Ti n = +/—, indicates spin up/down for the 3d 
electron of the Ti„ atom in the unit cell (see Fig. ED- 
Calculations were performed within the LCAO scheme 
with CRYSTAL09 and for the optimized structure. 


Magnetic configuration 

Etot 

E g 

a 

c 

AFM1 (+, -, +) 

0 

0.59 

5.167 

13.671 

AFM2 (+, -, +, -) 

6 




Nonmagnetic 

13 




FM 

13 




Expt.^ 


0.11 

5.157 

13.610 


a Experimental band gap from combined conductivity and 
thermoelectric coefficient measurements (RefJl3J). 
Experimental lattice parameters from Ref. Il35l . 
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TABLE IX: Reaction energies (in eV, estimated error 10 meV) for the TiC >2 reduction, computed with different 
hybrid functionals, and compared with experimental heats of reaction. The exchange fraction a = 1/coo (%) used in 
calculations on TiC >2 and Ti 2 03 is reported. The computed dielectric constants within PBEO and PBEOapBEo, 
respectively, are: 6.10 and 6.65 for TiC> 2 , and 5.63 and 7.81 for Ti 203 (estimated errors within 0.02). Calculations 

were performed within the LCAO scheme using CRYSTAL09. 




PBEO 

PBEOapBEo 

PBE0«pg EO 


Reaction 

a(Ti0 2 ) 

25 

16.4 

15.0 

Expt.^ 


a(Ti 2 Oa) 

25 

17.8 

12.8 


2 Ti 02 + H 2 —> Ti203 + H 2 O 


0.96 

1.22 

1.26 

1.30 

2Ti0 2 —> Ti 2 Q 3 + |0 2 


3.55 

3.81 

3.85 

3.81 


a Ref. 1138 and references therein. 


Ti 2 C >3 is a small gap semiconductor with a corundum 
structure. Its bulk unit cell is rhombohedral and con¬ 
tains two Ti 2 03 units, with atoms positioned as shown 
in Fig. [2] The magnetic properties of T^Os at low 
temperature have been subject of an intensive inves¬ 
tigation both at the experimenta l 137 : 138 and theoreti¬ 
cal leve h 55 : 56 ! 139 ! 140 However, no general consensus has 
been achieved on whether the low-temperature ground 
state of Ti 2 C >3 is antiferromagnetic (AFM) or nonmag¬ 
netic. We re-investigated the issue using both stan¬ 
dard and dielectric-dependent PBEO. Standard non¬ 
hybrid DFT predicts the material to be metallic^! at 
odds with the experimental evidence. Our spin-polarized 
PBEO calculations yielded an AFM ground state with 
AFM1 (+, —, —, +) spin configuration of Ti 3d elec¬ 
trons. Previous studies within Hartree-Fock^ hybrid 
DFT and DFT+U^ reported the same result. Ferro¬ 
magnetic (FM) and nonmagnetic solutions are found to 
be much higher in energy at the PBEO level. Inter¬ 
estingly, different investigations using screened-excliange 
hybrid functionals predicted the nonmagnetic configu¬ 
ration to be the most stable d 39 : 140 Calculations us¬ 
ing the dielectric-dependent PBEOapBEo functional con¬ 
firmed the magnetic configuration to be the AFM1, as re¬ 
ported in Table IVIIII The corresponding computed band 
gap of 0.59 eV is larger than the experimental value of 
0.11 eV obtained by combined conductivity and ther¬ 
moelectric coefficient measurementsbut still closer 
than our PBEO result of 2.63 eV, and in agreement 
with screened exchange hybrid calculations^^ reporting 



FIG. 2: Bulk primitive cell of corundum Ti 2 03 - Grey 
and red spheres represent Ti and O atoms, respectively. 


a value of 0.57 eV. 

For Ti02, we considered the rutile phase in the present 
analysis, since experimental data are available for the 
reaction enthalpy of the transformation to Ti 2 03 . We 
carried out full geometry optimizations within PBEO, 
PBEOapBEO) and partially self-consistent PBE0ap EEO . 
A band gap of 3.03 eV was obtained within PBE0apg EO , 
which is very close to our previously reported PBE0ap EE 
value of 2.96 eV (see Table IVTl) . confirming that the self- 
consistency is practically achieved. 

The computed energies of reactions © and © are 
summarized in Table IIXI 142 The dielectric-dependent 
PBEO yielded results in excellent agreement with ex¬ 
periment, with a clear improvement over standard 
PBEO. Previous first-principles investigations within 
GGA+ U 56 : 136 found that choosing the value of the Hub¬ 
bard parameter U between 2 and 3 eV results in reaction 
energies close to experiment. However, the chosen value 
of U strongly affects the computed electronic band struc¬ 
ture. Ti 2 03 was reported to be metallic at U = 2 eV, 
while U — 3 eV was needed to open a gap of compa¬ 
rable size to the experimental one. For TiC >2 it was 
demonstrated-^ that U should be set to ~ 10 eV in order 
to obtain a band gap of 3 eV. Another deficiency of the 
DFT+U approach lies in the failure to predict accurate 
lattice parameters for Ti 2 0 3 i££ Both of the above issues 
are elegantly solved within the dielectric-dependent hy¬ 
brid method, which we found out to be able to accurately 
characterize reaction energies, electronic structures and 
structural parameters within a single, parameter-free, 
and fully ab initio theoretical approach. 

IV. CONCLUSIONS 

We have reported a systematic analysis of band gaps, 
equilibrium geometries and phase stability of several 
ambient-pressure bulk polymorphs of wide gap oxide 
semiconductors of current interest for energy applica¬ 
tions. Hybrid xc functionals with a fraction of Hartree- 
Fock exchange weighted on the static electronic dielectric 
constant have been confirmed to provide an improvement 
over standard hybrid functionals concerning band gap es- 
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timation, the resulting accuracy being comparable to the 
more expensive GW method. Equilibrium geometries 
computed within the dielectric-dependent hybrid scheme 
have been found to be as accurate as with standard hy¬ 
brids. Predicted phase stability sequences have also been 
found to be the same with the two methods, with the ex¬ 
ception of WO3, in which a peculiar correlation between 
structure and band gap leads to different predicted stabil¬ 
ities. The stoichiometry transformation of TiC >2 to T^Oa 
has then been investigated. The dielectric-dependent hy¬ 
brid method has proved successful in describing band 
gaps and equilibrium geometries of both materials, while 
providing reaction energies for the Ti 02 reduction in ex¬ 
cellent agreement with experiments. In conclusion, the 
method qualifies itself for accurate determination of ei¬ 
ther excited- and ground-state properties of bulk solids, 


constituting a useful tool in those situations in which an 
accurate description of both is required within a common 
theoretical framework. 
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